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O 1 ABSTRACT 

O ' 



We use a systolic iV-body algorithm to evaluate the linear stability of the 
gravitational iV-body problem for iV up to 1.3 x 10 5 , two orders of magnitude 
greater than in previous experiments. For the first time, a clear iV-dependence of 
the perturbation growth rate is seen, fi e ~ IniV. The e-folding time for N = 10 5 



is roughly 1/20 of a crossing time. 
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Introduction 



Miller (1964) first noted the remarkable sensitivity of the gravitational iV-body prob- 
lem to small changes in the initial conditions. Errors or perturbations in the coordinates 
or velocities of one or more stars grow roughly exponentially, with an e-folding time that 
is of order the crossing time. The implication, verified in a number of subsequent studies 
(Lecar 1968; Hayli 1970), is that iV-body integrations are not reproducible over time scales 
that exceed a few crossing times. The instability is somewhat reduced when the Newto- 
nian force law is modified by a cutoff (Standish 1968), indicating that it is driven by close 

H | encounters. 
CZ i 

Of interest is the behavior of Miller's instability in the limit of large N. It is commonly 
assumed that the iV-body equations of motion go over to the collisionless Boltzmann equa- 
tion (CBE) as iV —>■ oo (Binney & Tremaine 1987). This would imply, for instance, that 
a particle trajectory which is integrable in a smooth potential should exhibit increasingly 
regular behavior as the number of point masses used to represent the smooth potential in- 
creases. On the other hand, if the rate of growth of small perturbations remains substantial 
even for large N, there would be an important sense in which the CBE does not correctly 
describe the behavior of iV-body systems. 

In fact there are indications that the growth rate of Miller's instability remains constant 
or even increases with N (Kandrup & Smith 1991; Heggie 1991; Goodman Heggie & Hut 1993), 
although this result is uncertain since published numerical experiments have been limited to 
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N < 10 3 . Here we describe the application of a new, "systolic" iV-body algorithm to this 
problem, which allows us to treat systems with N as large as 10 5 . We observe for the first 
time a clear iV-dependence of the instability: the growth rate is found to increase approxi- 
mately as IniV. Our methods and results are described in §2 and §3, and the implications 
for galactic dynamics are discussed in §4. 



2. Method 

Following Miller (1971), we integrated the coupled iV-body and variational equations: 

N 

Xj = —Gm ^ 

i=i 

N 

Xj = — Gm 

i=i 

Here Xj are the configuration-space coordinates of the ith particle and Xj are the components 
of its variational vector. The masses m are assumed equal. The variational equations 
represent the time development of the infinitesimal distance between two neighboring N- 
body systems. 

Equations (la) and (lb) were integrated using the systolic iV-body algorithm described 
by Dorband, Hemsendorf & Merritt (2002). This algorithm implements the fourth-order 
Hermite integration as described by Makino & Aarseth (1992). We adopted their formula 
for computing the time step of particle i, 



*(OI|a|%)| + |a?>fa)|*' 

Here a is the acceleration x, superscripts denote the order of the time derivative, t\ is the 
system time, and rj is a dimensionless constant; we set rj = 0.02. The same time step 
was used to integrate both the iV-body and variational equations. The systolic algorithm 
distributes the particles equally among p processors and computes forces by systematically 
shifting the particle coordinates between processors in a ring. A single processor was used 
for small particle numbers while 64 processors were used for the largest-iV integrations 
(Table 1). The multi-processor integrations were carried out using the Cray T3E at the 
Hochstleistungsrechenzentrum in Stuttgart. 

Initial conditions were generated randomly from the isotropic Plummer model, whose 
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Fig. 1. — Growth of the variation in the spatial coordinate for 10 Plummer-model integrations 
with iV = 32768. (a) A g , the geometric mean of the variations, (b) A a , the arithmetic mean 
of the variations, (c) A m , the maximum variation. 
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density and potential satisfy 

, , 3GM b 2 GM 
p( r ) = — — , $ r = 3 

We adopted standard iV-body units (Heggie & Mathieu 1986) such that G = M = 1, E = 
— 1/4, giving a scale factor b = 3n/lQ. We defined the crossing time t cr as R/V, with 
R = -GM 2 /2E and V 2 = -2E/M; in these units, t cr = 2-\/2. The variational vectors X 
and V were assigned an initial amplitude of 10~ 30 for each particle with randomly chosen 
directions. 

The parameters of the integrations are listed in Table I. Each integration was continued 
until a time of 20, or roughly 7 crossing times, with the exception of the largest run, N = 
131072 which was terminated at t — 12. This time is short enough that two-body relaxation 
should not be important except perhaps for the smallest N, and long enough to show a 
clearly exponential growth of the solutions of the variational equations. 

We focussed on the amplitudes of the Cartesian components of the variational vectors, 
Xj = (Xj, Yi, Zi), since the variational velocities tend to exhibit spikes in their time depen- 
dence (Miller 1964). We examined three choices for the amplitude A of the separation: 

1. A m , the maximum over % of the |Xj|. 

2. A a , the arithmetic mean of the |Xj|. 

3. A g , the geometric mean of the |Xj|. 

The instability growth rate \i e was defined as 

p. = "^fa'-'"^''. (4) 
£2 — ti 

Except in the case of small N, N < 1024, growth in the variation was found to be nearly 
exponential and hence the computed values of fj, e depended only weakly on t 2 and t 1 . We 
chose ti — 1 and t 2 = 20, except for the N = 131072 run for which t 2 = 12. 



3. Results 

The three amplitudes A defined above were found to be very similar in most of the 
integrations, as shown in Figure 1 for the 10 integrations with N = 32768. In what follows 
we adopt A = A g , the geometric mean. 

Figure 2a shows t e = /i" 1 for each of the integrations. The mean value of the e-folding 
time and its uncertainty are plotted in Figure 2b; the latter was defined as the standard 
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Fig. 2. — Exponentiation times for the iV-body integrations, (a) t e for each of the integra- 
tions. The open circle is the mean value of t e quoted by Goodman et al. (1993) from seven 
integrations of an iV=256 Plummer model, (b) Mean values of t e for each N. Error bars 
represent the standard error of the mean. 
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error of the mean, or n~ l l 2 times the standard deviation, with n the number of distinct 
iV-body integrations. For the first time, a clear iV-dependence can be seen, in the sense that 
the average value of t e declines with increasing N: large- N systems are more unstable than 
small- iV systems. 

A number of predictions have been made for the large- N dependence of the e-folding 
time. Gurzadyan & Savvidy (1986) estimated t e ~ N l ^t cr based on a geometrical approach. 
This iV dependence is clearly inconsistent with Figure 2. Gurzadyan & Savvidy 's approach 
was criticized already by Heggie (1991) and Goodman, Heggie & Hut (1993) due to its 
improper treatment of close encounters. The latter authors argued that t e /t cr ~ 1/lniV or 
~ l/ln(lniV); the weaker dependence would hold only after a time long enough that the 
perturbation from one star was able to propagate through the system. 

We tested these predictions against the iV-body data. Figure 3 shows fits of two func- 
tional forms to the mean growth rates: 



Since any such relation is expected to be valid only in the limit of large N, we restricted 
the fits to iV > 1024. The best-fit parameters are given in Table 2. We used a standard 
least-squares routine that accounts for errors in the dependent variable (/x e ); in the case of 
the data pointa with iV = 65536 and 131072, for which there was only one integration, the 
uncertainty in \i e was assumed to be the same as in the integration with iV = 32768. These 
two data points were omitted when computing the values of x 2 given in Table 2. 

While /jL e t cr oc iV° dependence can clearly be ruled out, we can not distinguish between 
a IniV and a ln(lniV) dependence. Expressed in terms of t e and t cr , the best-fit relations 
are 



He = a + b In N, 
fi e = c + d\n(\nN). 



(5a) 
(5b) 



te/t, 



0.99 



(6a) 



cr 



In (1.09 x 10 3 iV)' 
0.112 



te/t, 



cr 



In (0.70 In TV)' 



(6b) 



Distinguishing between these two functional forms, or other similar ones, would clearly be 
very difficult; even for iV = 10 6 , the two relations predict values of t e that differ only by 
-5%. 





Fig. 3. — Fits of the iV-body growth rates to two functional forms, /i e oc IniV and /i e oc 
ln(lniV). Fitting parameters are given in Table 2. 
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4. Discussion 

Our results confirm that the gravitational TV-body problem is inherently chaotic and 
furthermore that the degree of chaos, as measured by the rate of divergence of nearby 
trajectories, increases with increasing N. We have directly established a characteristic e- 
folding time of ~ t cr /20 for systems with N 10 5 ; if the weak iV-dependence found here 
can be extrapolated, the characteristic time for systems containing ~ 10 12 particles would 
be only slightly smaller, ~ t cr /30. Hence trajectories in stellar and galactic systems diverge 
on a time scale that is generically much shorter than the crossing time. 

If the collisionless Boltzmann equation (CBE) is a valid representation of large- N stellar 
systems, it should be possible to show that the iV-body trajectories go over, in the limit of 
large N, to the orbits in the corresponding smoothed-out potential. The generic instability of 
the iV-body problem precludes this, since the characteristics of the CBE can not be identified 
with the iV-body orbits for times longer than ~ t cr . At the same time, much experience 
with iV-body integrations demonstrates that in many ways the behavior of large- N systems 
matches expectations derived from the CBE (e.g. Aarseth & Lecar 1975). 

A likely resolution of this seeming paradox is that the macroscopic, or finite-amplitude, 
behavior of trajectories is not well predicted by the rate of growth of small perturbations. 
For instance, orbits integrated in "frozen" iV-body potentials behave more and more like 
their smooth-potential counterparts as iV is increased, even though their Liapunov expo- 
nents remain large (Valluri & Merritt 2000; Kandrup & Sideris 2001). The initial growth of 
perturbations is exponential but it saturates, at an amplitude that varies inversely with N 
(Valluri & Merritt 2000; Hut & Heggie 2001). Thus the CBE may be a good predictor of 
the macroscopic dynamics of large- N systems even if it does not reproduce the small-scale 
chaos inherently associated with iV-body dynamics. 

This work was supported by NSF grant 00-71099 and by NASA grants NAG5-6037 
and NAG5-9046 to DM. We thank S. Aarseth, J. Goodman, D. Heggie and H. Kandrup 
for useful discussions. We are grateful to the Hochstleistungsrechenzentrum in Stuttgart for 
their generous allocation of computer time. 
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Table 1. Parameters of the iV-body integrations. 
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128 


175 


1 


0.243 


0.005 


4.41 


0.10 


256 


175 


1 


0.212 


0.004 


4.93 


0.08 


512 


175 


1 


0.209 


0.002 


4.89 


0.06 


1024 


175 


1 


0.201 


0.002 


5.04 


0.04 


2048 


175 


1 


0.195 


0.005 


5.26 


0.04 


4096 


35 


1 


0.185 


0.003 


5.43 


0.078 


8192 


10 


1 


0.171 


0.002 


5.86 


0.078 


16384 


10 


64 


0.169 


0.001 


5.92 


0.038 


32768 


10 


64 


0.160 


0.002 


6.28 


0.086 


65536 


1 


64 


0.156 




6.42 




131072 


1 


64 


0.143 




6.99 





Note. - N is the particle number; n is the 
number of distinct integrations; p is the number 
of processors used; (t e ) and a te are the mean e- 
folding time and its error; (/i e ) and a^ e are the 
mean e-folding rate and its error. 



Table 2. Results of the linear regression fits, Y = a + bX. 



X, Y Variables 


a 


P 




In N, n e 


2.51 ±0.12 


0.36 ±0.01 


2.0 


In (In iV), n e 


-1.11 ±0.26 


3.14 ±0.12 


3.3 



